%% 4th order runge-kutta method (integration) to solve the 
%% ordinary differential equation (ODE) approximative:
% input:
%   fh       ...  function handle
%   t0       ...  initial value of the independent variable.
%   tf       ...  final value of the independent variable.
%   y0       ...  initial value(s) of the dependent variable(s).
%                 to solve a system of equations, all initial values should
%                 be contained in a row vector.
%   iNSteps  ...  number of time steps between "t0" and "tf".
%
% output:
%   ta       ... vector containing all values of the independent variable
%                at which the approximated solutions has been obtained.
%   ya       ... vector/matrix containing approximated solution values of
%                the ODE.

%function [ta, ya] = solveODE_RK4(fh, t0, tf, y0, ki, iNSteps)
function [ta, ya] = solveODE_RK4(fh, t0, tf, y0, iNSteps)
    % initialize ...
    iLen = iNSteps + 1;
    iNEqs = length(y0); % number of equations ...
    ta = linspace(t0, tf, iLen);
    ya = zeros(iNEqs, iLen);
    
    % calculating the step size ...
    h = abs(tf - t0)/iNSteps;
    
    % given initial value problem: y'(x) = f(x, y(x)) and
    % initial condition:
    ya(1:iNEqs, 1) = y0';
    %ta(1) = t0;
    
    for i = 1:iNSteps
        % calculating the slopes (slope = weighted average of slopes) ... 
%         k1 = feval( fh, t0, y0, ki );
%         k2 = feval( fh, t0 + h/2.0, y0 + (h/2.0)*k1, ki );
%         k3 = feval( fh, t0 + h/2.0, y0 + (h/2.0)*k2, ki );
%         k4 = feval( fh, t0 + h, y0 + h*k3, ki );

        k1 = feval( fh, t0, y0 );
        k2 = feval( fh, t0 + h/2.0, y0 + (h/2.0)*k1 );
        k3 = feval( fh, t0 + h/2.0, y0 + (h/2.0)*k2 );
        k4 = feval( fh, t0 + h, y0 + h*k3 );
        % applying the 4th order Runga-Kutta formula ...
        y0 = y0 + (k1 + 2.0*k2 + 2.0*k3 + k4)*h/6.0;
        % adding the step size ...
        t0 = t0 + h ;
        
        % save the result(s) ...
        ya(1:iNEqs, i+1) = y0';
    end    
end
